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We develop an extension of the Gutzwiller approximation to finite temperatures based on the 
Dirac-Frenkel variational principle. Our method does not rely on any entropy inequality, and is 
substantially more accurate than the approaches proposed in previous works. We apply our theory to 
the single-band Hubbard model at different fillings, and show that our results compare quantitatively 
well with dynamical mean field theory in the metallic phase. We discuss potential applications of 
our technique within the framework of first principle calculations. 
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The Gutzwiller approximation (GA)P10 is a very use¬ 
ful tool in order to study the ground state of complex 
strongly correlated electron systems. This important 
many-body technique has been also formulated and im¬ 
plemented in combination with density functional the¬ 
ory (DFT),® e.g., in the LDA-I-GA approach,ISlt^ which 
has been applied successfully to many real materials 
For strongly correlated metals, the accuracy of the 
GA is c ompa rable with dynamical mean field theory 
(DMFT), b^F^ I even though the GA is much less compu¬ 
tationally demanding. This property makes it an ideal 
theoretical tool, as numerical speed is essential for the 
purpose of studying and discovering new materials. 

In order to study several temperature-dependent phe¬ 
nomena, such as structural and magnetic transitions 
and coherence-incoherence crossovers, it would be highly 
desiderable to have at our disposal an extension to finite 
temperatures of the GA as accurate as the ordinary the¬ 
ory for the ground state. In fact, this would enable us 
to study these properties also for correlated systems so 
complex to be out of the reach of the presently available 
methods, such as DMFT. 

An extension of the GA to finite temperatures has 
been previously proposed in Refs. 1181191 This approx¬ 
imation scheme is based on an exact entropy inequal¬ 
ity which enables to calculate an upper bound to the 
free energy,li^ and minimize it numerically. Of course, 
underestimating the entropy using an entropy inequal¬ 
ity — rather than calculating it exactly — constitutes 
a source of approximation not present in the ordinary 
zero-temperature GA. In particular, it has been shown 
that this additional source of approximation generates a 
few pathologies of the theor y, suc h as giving a negative 
entropy at low temperatures.^i^*^ 

In this work we introduce an extension of the GA to fi¬ 
nite temperatures based on the Dirac-Frenkel variational 
principle^ ^"^ a nd, in particular, on the time-dependent 
GA theorji^Hm (that we generalize to mixed states). Our 
method does not rely on any entropy inequality, but only 
on the variational principle and the Gutzwiller approx¬ 
imation — which are the same approximations done in 
the ordinary zero-temperature GA. Consequently, as we 


are going to show, our theory improves considerably the 
method of Refs. mm and gives results in good quanti¬ 
tative agreement with DMFT for correlated metals, even 
though it is much less computationally demanding. 

Imaginary-time evolution. — Let us consider a generic 
system of correlated electrons represented by a Hamil¬ 
tonian Pi, and dehne the imaginary-time evolution of a 
given initial density matrix po as follows: 

p(r)=e-«"poe-«V (1) 

i.e., according to the following differential equation: 

9rP(T) = -{Hpir) + p{T)n) = -{PL, /3(r)} . (2) 

Our aim consists in approximating the imaginary-time 
dynamics defined above and use it to construct the state 
of N electrons at temperature T. In fact, if r = /3/2 
and po = Pn is the projector onto the subspace with N 
electrons, Eq. ([^ reduces to which represents 

a thermal state with T = 1/ 

In order to derive our approximation scheme, it will be 
useful to think of p as the density matrix corresponding 
to an ensemble of pure states 

n 

where are fixed probabilities coefficients. Within this 
definition, evolving p according to Eq. Q amounts to 
evolve all of the pure states of the ensemble according to 
the equation 

d|^'„(r)) = -?i|4'„(r)) dr . (4) 

Note that Eq. Q resembles a Schrodinger evolution in 
imaginary time, as it can be obtained from the ordinary 
real-time Schrodinger evolution 

d\'^n(t)) = -iPi\A>n{t)) dt (5) 

by substituting dt —> —i dr. 

Real-time Dirac-Frenkel scheme .— Let us introduce 
the following action!^ 

5{p„}[{4'„(<)}] = ^ dtC[p^}[{'^n{t)}] (6) 
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n 

which depends parametrically on the probability coef¬ 
ficients Pn (that are fixed). From now on we refer to 
Eq. Q as the Dirac-Frenkel action. It can be readily 
verified that, regardless the values of the exact solu¬ 
tion of the Lagrange equations for the ensemble of states 
{|5'„(t))} is given by Eq. ([5|. 

The key advantage of the Dirac-Frenkel characteriza¬ 
tion of the time evolution outlined above is that it allows 
us to build up a well-founded variational approximation 
scheme for the time evolution [Eq. ([^] as follows. 

Let us assume that we want to solve approximately 
the time-dependent problem by restricting the search of 
the solution within a submanifold TW of trial ensembles 
{IOnce we are able to evaluate the action S along 
any given trajectory in A4, the Dirac-Frenkel variational 
principle provides us with a prescription to approximate 
the instantaneous time evolution of any {['1'™)} € -M. 
Note that, by construction, this time evolution is such 
that {l'I'n(t))} G M. \/t. 

Application to the GA .— For sake of simplicity, in this 
work the method will be formulated for the single-band 
Hubbard model: 

k cr=t,t R 

where k is the momentum conjugate to the site label 
R and a is the spin label. The extension to multi-band 
Hubbard models is straightforward, and its numerical im¬ 
plementation will be discussed in a future work. In order 
to benchmark our theory, we present finite-temperature 
calculations of the Hamiltonian [Eq. ([^] at different fill¬ 
ings N/Af = 1 -I- (5, where JV is the number of /c-points 
and S is the doping. 

Here we want to search for the saddle point of the 
Dirac-Frenkel action within the manifold of ensem¬ 
bles of Gutzwiller states represented as follows: 

mn)} = {rG\^0n)}=MG, (9) 

where |5'o„) are Slater determinants and Vg = YIr'^r 
is an operator whose local components are defined as 
Vr = Ar |i?, F)(i?, r|, where Ar are numbers and 


|i?, F) (i?, r| are the projectors onto the corresponding lo¬ 
cal many-body states |i?,F) S {[0), ti)}- 

The physical density matrix corresponding to the en¬ 
semble [Eq. ([^] is pG = Vq Po Vq, where 

PI = |«'0n)(^0rx| / ^ Pn (4'Onl^'On) (10) 

n n 

is called variational density matrix. We assume that pj 
can be represented as the Boltzmann distribution of a 
generic noninteracting Hamiltonian V t. In order to cal¬ 
culate the energy corresponding to pG — which is nec¬ 
essary to evaluate the Dirac-Frenkel action, see Eq. Q, 
— the manifold of ensembles A4 g is further restricted by 
the so called Gutzwiller constraintsP^^l^IlII 

Tr[pSPl^«] = l (11) 

Tr[pS K'Pr cLcfl.] = Tr[pS = [1 + 6]/2 . (12) 

Furthermore, the GA is assumed, which is an approxi¬ 
mation scheme that, as DMFTjC^ becomes exact in the 
limit of infinite coordination lattices. 

As in Ref. nzi we introduce the matrix of slave-boson 
amplitudes: 

(^PP, = (5pp/Ar (13) 

PpO = Tr[p5|i?,r)(R,F|] . (14) 

Within the above def inition s, the Gutzwiller constraints 
can be represented as! ^^ * ^^ l 

Tr[0l^] = l (15) 

Tr[0t0FtFj = Tr[pS cLc^.] = [1 + S]/2 , (16) 

where [iG]rr' = (r 1"^) ■ Furthermore, it can be 
shown that represents the local reduced density ma¬ 
trix in the basis {|i?, F)}, while the expectation values of 
quadratic non-local observables is given by: 

Tr[PG cLcflv] = Tr[p* c^c^, J , (17) 

where TZ = Ti[(j)^F^(j>F^]/[l — Using the above 

equations, the GA Dirac-Frenkel Lagrange function can 
be rewritten as follows!^ 


n k cr=td 

+TT[(l>'<idt(l>] -Tr U(j)cj)'< F^F^fIf^ - ^ (Ti[D (/)"< F^^(I)F^] - VTl [l - + c.c.'^ . 

o-=td 

I 


Note that, following Ref. [3 we have formally enforced the definition of TZ using the Lagrange multiplier T>. 
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The Lagrange equations for the real-time dynamics in¬ 
duced by Eq. (19) are the following: 


i^-on) = o Vn 

[idt-Hf^^[V,V*]] </, = 0 

7^ = Tr[<^t^t^^^] 


V = 2[l-5^ 


Tr 




where 


^ |7^P5]efecLc,, 


ka 


Hf^^[D,V*]cl> = ^{Tr 


U( 


At pt 


E'KF F, 


+ ^ (Tr[l?</.tFt<^E^] +c.c.)}</> 


(19) 

( 20 ) 
( 21 ) 
( 22 ) 

(23) 

(24) 
■ (25) 


Note that the generator of the instantaneous evolution is 
quadratic and identical for all of the |'I'o„), and that also 
the evolution of (j) resembles formally a time-dependent 
Schrodinger equation. 

The instantaneous real-time evolution described by the 
equations above corresponds to apply well defined incre¬ 
ments on all of the the states of AIg, see Eq. (|^. We 
may represent these increments as follows: 

d|^„) = [{dtVG) l^-On) + Vg (Stl^'On))] ctt . ( 26 ) 

Imaginary-time dynamics. — Our goal consists in mod¬ 
ifying the real-time GA dynamics defined above in order 
to approximate the imaginary-time evolution [Eq. Q] . 

The formal similarity between Eqs. Q and ([^ suggests 
us that it is possible to approximate the imaginary-time 
evolution of simply by substituting dt —A —idr 


in Eq. (261. It can be readily verified that this prescrip¬ 
tion would amount to update the Gutzwiller variational 
parameters, see Eqs. (13) and (14), as follows!^ 


dr + n^p[n,n*]\ j^-on) = o v 

[dr + Hf^^[v,v*]] = 0. 


(27) 

(28) 


Unfortunately, Eqs. (27) and (28) violate the Gutzwiller 
constraints, see Eqs. (15) and (16). Gonsequently, sim¬ 
ilarly to Ref. HBl it is necessary to define a “projection 
scheme” in order to enforce them at every time step. 

Here we propose to enforce Eqs. (15) and (16) by using 
the following prescription: 

d^+H%[TZ,n*,Eo]\ l^-on) = 0 Vn (29) 

[dr + Felb P. F*, A^ E% = 0, (30) 

where the “generators” have been modified as follows: 
nZ = - Eo (31) 


H. 


' = Ef^^4>+ 


embV 


54>^ 


(^,(32) 


and Eq{t) is constructed in order to enforce the normal¬ 
ization condition of Pq, see Eq. (10), while E‘^{t) and 
A‘^(t) are constructed in order to enforce Eqs ( flSj ) and 
(16), respectively. 


We point out that the procedure defined above enables 
us to recover the ordinary GA theory for the ground state 
at T —)■ oo. In fact, within the formulation of Ref. [71 the 
GA parameters of the ground-state are obtained as the 
ground states of 77qp and 77^^, which correspond to a 
fix point of our imaginary-time dynamics. 

It can be readily verified that Eq. (29) implies that 
the imaginary-time evolution of the variational density 
matrix is given by: 

P*{t) = ^ ^ 33 ^ 

where Z{t') = \TZ{t')\‘^ is the Gutzwiller quasi-particle 
weight, and Fo(t') is constructed in order to enforce the 
normalization condition of Pq{t) for all imaginary times. 
In fact, Eq. (33) satisfies: 


9tpUt) = / 5 o(t)} , 


(34) 


which is consistent with Eq. ( |30| ), and enables us to avoid 
to keep track of the time evolution of all of the states of 
Mg (which would be practically impossible). 

Note that, since we are in the thermodynamical limit, 
the expectation values with respect to Pq{t) can be evalu¬ 
ated in the grand-canonical ensemble, i.e., we can assume 
that 


PoM 




(35) 


where /35 (t) = 2 JJ dT'Z(T'), N is the number operator, 
and Pq{t) is such that the system has N electrons in 
average. 

The imaginary-time evolution of the slave-boson am¬ 
plitudes is obtained by substituting Eq. (351 into the La¬ 


grange equations for 0, A'^, TZ, V, F° and solving them nu¬ 
merically. 

Numerical results. — Let us now discuss our numerical 
calculations of the Hubbard model, see Eq. (|^. We as¬ 
sume a semicircular density of states (corresponding to 
a Bethe lattice in infinite dimensions)^ and set the half¬ 
bandwidth D as the unit of energy. Eor comparison, we 
perform DMFT calculations using the continuous time 
quantum Monte Carlo method with hybridization expan- 
sioipnias impurity solver, as implemented in TRIQS.^^ 

In the upper panel of Fig. f^is shown the evolution of 
the double occupancy d = as a function 

of the temperature at half-filling for several values of U. 
In the lower panel is shown the corresponding evolution 
of the total energy £. The GA results are shown in com¬ 
parison with DMFT and the Gutzwiller data of Ref. HHI 

The agreement between the GA and DMFT-I-CTQMC 
is quantitatively satisfying, especially for smaller values 
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Ar=i (5=0) 



FIG. 1: GA calculations of the single-band Hubbard model 
at half-filling {N = 1) in comparison with DMFT-I-GTQMC 
and the data of Ref. 1181 Upper panel: evolution of the double 
occupancy as a function of the temperature. Lower panel: 
evolution of the total energy as a function of the temperature. 


of U and higher temperatures (i.e., when the system 
is less correlated). Indeed, our method improves sub¬ 
stantially the results obtained within the approximation 
scheme of Ref. m The slight quantitative discrepancy for 
larger C/’s reflects the known fact that the Mott insulator 
is not well described by the GA, but is approximated by 
the simple atomic limit — that is a state with d = 0. 
However, as long as the system is metallic, our extension 
of the GA to finite temperatures is remarkably accurate. 

Let us now consider the Hubbard model away from 
half-filling. In particular, we consider the case of iV = 0.8 
electrons per site (i.e., 6 = —0.2). In the upper panel of 
Fig .[^is shown the temperature dependence of the double 
occupancy for several values of U, while in the lower panel 
is shown the evolution of the total energy £. Finally, in 
the inset of the lower panel is shown the temperature 
dependence of the entropy for t//Z/ = 4, in comparison 
with the DMFT entropy calculated in Ref. [29l 

We point out that, as discussed before, the entropy is 
not evaluated directly from the GA variational param¬ 
eters (which could be done only approximately, e.g., by 


Ar=0.8 ((5 = -0.2) 



FIG. 2: GA calculations of the single-band Hubbard model 
away from half-filling {N = 0.8) in comparison with 

DMFT-I-GTQMG. Upper panel: temperature dependence of 
the double occupancy. Lower panel: temperature dependence 
of the total energy £{T). Inset of the lower panel: tempera¬ 
ture dependence (in logaritmic scale) of the GA entropy S in 
comparison with the DMFT data of Ref. 1291 


using the entropy inequality of Ref. [T9|), but is calcu¬ 
lated from the imaginary-time evolution of the total en¬ 
ergy using the well known thermodynamical identities 
dS = d£/T, S{T = 0) = 0. Note that the value of S at 
T —>• oo calculated according to these equations depends 
on the evolution of the total energy within the whole 
range of temperatures. It is for this reason that the GA 
entropy shown in Fig. is slightly shifted with respect to 
DMFT at high temperatures — even though the atomic 
limit belongs to the GA variational space, and is thus 
captured exactly by our approximation scheme.l^ 

The agreement between the GA and DMFT-I-CTQMG 
is even better for N = 0.8 than for half-filling (which is to 
be expected, as the doped system is metallic for all C/’s). 
In particular, it is remarkable that the agreement for S 
is satisfying for U/D = 4, which is the largest interaction 
strength considered. 

In conclusion, we have developed an extension of the 
Gutzwiller approximation to finite temperatures based 
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on the Dirac-Frenkel variational principle. Since our 
method does not rely on any entropy inequality, but 
only on the variational principle and the Gutzwiller ap¬ 
proximation, it is as accurate as the ordinary GA the¬ 
ory for the ground state, and improves substantially the 
method previously proposed in Refs. 11811^ We have per¬ 
formed benchmark calculations of the single-band Hub¬ 
bard model at different hllings, and compared our results 
with DMFT-I-GTQMG, finding good quantitative agree¬ 
ment between the two methods in the metallic phase. We 
believe that our method will enable us to calculate from 
first principles several important physical quantities — 
such as the specific heat, the entropy and the tempera¬ 


ture dependent structural properties — of strongly cor¬ 
related systems presently too complex to be studied with 
more accurate methods, such as DMFT. 
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Note that DMFT is an exact theory for this system. 

The reason why DMFT does not suffer this inconvenience 
is that it is an exact theory in infinite dimensions (while 
the GA is a variational approximation). 





